namespace NBody.Analysis
{
    using System;
    using NBody.Cosmology;
    
    public class VirialRadius
    {
        public static double CalculateR200(NBodySystem s, double background_density)
        {            
            double rho200 = 200.0 * background_density;
            Console.Error.WriteLine("Finding virial radius for rho_200 = {0}", rho200);

            if (!s.IsSortedByRadius)
                s.Sort();
            
            int i = 0;
            double m_tot = s[i].Mass;
            double r = s[i].Radius;
            double rho = m_tot / (4.0 / 3.0 * Math.PI * r * r * r);
            
            while (rho > rho200)
            {
                i++;
                m_tot += s[i].Mass;
                r = s[i].Radius;
                rho = m_tot / (4.0 / 3.0 * Math.PI * r * r * r);
            }
            
            Console.Error.WriteLine("  rho   = {0}", rho);
            Console.Error.WriteLine("  r_200 = {0}", r);
            Console.Error.WriteLine("  M_200 = {0}", m_tot);
            Console.Error.WriteLine("  N_200 = {0}", i);
            
            Console.Error.WriteLine("  Powers et al. suggest a softening of {0}", 4.0 * r / Math.Sqrt((double)i));
            
            return r;
        }
        
        public static double CalculateR200(NBodySystem s, Model model)
        {            
            double rho200 = 200.0 * model.Density(model.Redshift(s.Time));
            Console.Error.WriteLine("Finding virial radius for rho_200 = {0}", rho200);

            if (!s.IsSortedByRadius)
                s.Sort();
            
            int i = 0;
            double m_tot = s[i].Mass;
            double r = s[i].Radius;
            double rho = m_tot / (4.0 / 3.0 * Math.PI * r * r * r);
            
            while (rho > rho200)
            {
                i++;
                m_tot += s[i].Mass;
                r = s[i].Radius;
                rho = m_tot / (4.0 / 3.0 * Math.PI * r * r * r);
            }
            
            Console.Error.WriteLine("  rho   = {0}", rho);
            Console.Error.WriteLine("  r_200 = {0}", r);
            Console.Error.WriteLine("  M_200 = {0}", m_tot);
            Console.Error.WriteLine("  N_200 = {0}", i);
            
            Console.Error.WriteLine("  Powers et al. suggest a softening of {0}", 4.0 * r / Math.Sqrt((double)i));
            
            return r;
        }
    }
}
            
